!***********************************************************************************************************8
! This chi object provides a fit to the distance function chi(z)
!
! Declare:
!
! type(tChi) :: chi
!
! call new(chi, [universe])
!
! with optional user-supplied universe type.
!
! Now you can use a_chi(chi, a)
!
! to calculate the distance to scale factor a = 1 / (1 + z)
!***********************************************************************************************************8

module mChi

        use mPrecision
        use mObject
        use mUniverse
        use mFit

        implicit none

        real(D), parameter :: LL_A = 0.0009_D, &
                              UL_A = 1.0_D

        interface new
                module procedure new_chi
        end interface new

        interface getUniverse
                module procedure getUniverse_chi
        end interface getUniverse

        interface a_chi
                module procedure a_chi_chi
        end interface a_chi

        interface getH
                module procedure getH_chi
        end interface getH

        interface getOmega_m0
                module procedure getOmega_m0_chi
        end interface getOmega_m0

        interface setOmega_b
                module procedure setOmega_b_chi
        end interface setOmega_b

        interface setOmega_d
                module procedure setOmega_d_chi
        end interface setOmega_d

        interface setOmega_de
                module procedure setOmega_de_chi
        end interface setOmega_de

        interface setW_0
                module procedure setW_0_chi
        end interface setW_0

        interface setW_a
                module procedure setW_a_chi
        end interface setW_a

        interface evaluate
                module procedure evaluate_chi
        end interface evaluate

        private
        public :: new
        public :: getUniverse
        public :: getH
        public :: getOmega_m0
        public :: setW_0, setW_a
        public :: a_chi
        public :: chiWriteCoeffs
        public :: LL_A, UL_A

contains

        subroutine new_chi(c, universe)
                type(tChi), intent(out) :: c
                type(tUniverse), intent(in), optional :: universe
                if( present(universe) ) then
                        c%u = universe
                else
                        call new(c%u)
                end if
                call new(c%chiFit, LL_A, UL_A, universe = c%u)
                call fit(c%chiFit, chi)
        end subroutine new_chi

        real(D) function chi(fitParams, a)
                use mConstants
                use mIntegration
                use mFit
                type(tFitParams), intent(in) :: fitParams
                real(D), intent(in) :: a
                type(tIntegral) :: chiIntegral

                call new( chiIntegral, a, 1.0_D, universe = getUniverse(fitParams) )
                call integrate(chiIntegral, chiIntegrand, chi)
!                if( .not. converged(chiIntegral) ) &
!                        write(0, *) 'WARNING: chi failed to converge'
                call checkConvergence(chiIntegral, 'chi failed to converge.')
                chi = c * chi
        end function chi

        real(D) pure function chiIntegrand(integral, x)
                use mIntegration
                type(tIntegral), intent(in) :: integral
                real(D), intent(in) :: x ! I changed variables x = 1 / a
                type(tUniverse) :: u
                u = getUniverse(integral)
                chiIntegrand = 1.0_D / ( hubble(u, x) * x ** 2 )
        end function chiIntegrand

        subroutine chiWriteCoeffs(c)
                type(tChi), intent(in) :: c
                call writeRFitCoeffs(c%chiFit)
        end subroutine chiWriteCoeffs
        
        type(tUniverse) pure function getUniverse_chi(c)
                type(tChi), intent(in) :: c
                getUniverse_chi = c%u
        end function getUniverse_chi

        real(D) function getH_chi(c)
                type(tChi), intent(in) :: c
                getH_chi = getH( getUniverse(c) )
        end function getH_chi

        real(D) function getOmega_m0_chi(c)
                type(tChi), intent(in) :: c
                getOmega_m0_chi = getOmega_m0( getUniverse(c) )
        end function getOmega_m0_chi

        real(D) pure function ll(c)
                type(tChi), intent(in) :: c
                ll = fitLL(c%chiFit)
        end function ll

        real(D) pure function ul(c)
                type(tChi), intent(in) :: c
                ul = fitUL(c%chiFit)
        end function ul

        real(D) function chiLL(c)
                type(tChi), intent(in) :: c
                chiLL = a_chi(c, ll(c))
        end function chiLL

        real(D) function chiUL(c)
                type(tChi), intent(in) :: c
                chiUL = a_chi(c, ul(c))
        end function chiUL

        real(D) pure function evaluate_chi(chi,x)
                type(tChi), intent(in) :: chi
                real(D), intent(in) :: x 
                evaluate_chi = evaluate(chi%chiFit, x)
        end function evaluate_chi

        real(D) function a_chi_chi(c, a)
                type(tChi), intent(in) :: c
                real(D), intent(in) :: a
                if(a < LL_A .or. a > UL_A) then
                        write(0, *) 'ERROR: a_chi called with a out of range:'
                        write(0, *) '     ', LL_A, '<=', a, '<=', UL_A
                        stop
                end if
                a_chi_chi = evaluate(c%chiFit, a)
        end function a_chi_chi

        subroutine chiCheckU(c, u)
                type(tChi), intent(in) :: c
                type(tUniverse), intent(in) :: u
                if( .not. equal(u, c%u) ) then
                        write(0, *) 'ERROR: wrong chi fit.'
                        stop
                end if
        end subroutine chiCheckU

        subroutine setOmega_b_chi(chi, omega_b)

                type(tChi), intent(inout) :: chi
                real(D), intent(in) :: omega_b

                type(tUniverse) :: u

                u = getUniverse(chi)
                call setOmega_b(u, omega_b)
                call new(chi, universe = u)

        end subroutine setOmega_b_chi

        subroutine setOmega_d_chi(chi, omega_d)

                type(tChi), intent(inout) :: chi
                real(D), intent(in) :: omega_d

                type(tUniverse) :: u

                u = getUniverse(chi)
                call setOmega_d(u, omega_d)
                call new(chi, universe = u)

        end subroutine setOmega_d_chi

        subroutine setOmega_de_chi(chi, Omega_de)

                type(tChi), intent(inout) :: chi
                real(D), intent(in) :: Omega_de

                type(tUniverse) :: u

                u = getUniverse(chi)
                call setOmega_de(u, Omega_de)
                call new(chi, universe = u)

        end subroutine setOmega_de_chi

        subroutine setW_0_chi(chi, w_0)

                type(tChi), intent(inout) :: chi
                real(D), intent(in) :: w_0

                type(tUniverse) :: u

                u = getUniverse(chi)
                call setW_0(u, w_0)
                call new(chi, universe = u)
                
        end subroutine setW_0_chi

        subroutine setW_a_chi(chi, w_a)

                type(tChi), intent(inout) :: chi
                real(D), intent(in) :: w_a

                type(tUniverse) :: u

                u = getUniverse(chi)
                call setW_a(u, w_a)
                call new(chi, universe = u)
                
        end subroutine setW_a_chi

end module mChi
